Reservoir engineering and dynamical phase transitions in optomechanical arrays 
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We study the driven-dissipative dynamics of photons interacting with an array of micromechanical 
membranes in an optical cavity. Periodic membrane driving and phonon creation result in an 
effective photon-number conserving non-unitary dynamics, which features a steady state with long- 
range photonic coherence. If the leakage of photons out of the cavity is counteracted by incoherent 
driving of the photonic modes, we show that the system undergoes a dynamical phase transition 
to the state with long-range coherence. A minimal system, composed of two micromechanical 
membranes in a cavity, is studied in detail, and it is shown to be a realistic setup where the key 
processes of the driven-dissipative dynamics can be seen. 



PACS numbers: 42.50.Wk, 64.70.Tg, 42.50.Lc 

I. INTRODUCTION 

While quantum computing and quantum simulation 
are traditionally discussed as dynamics of isolated many- 
body systems governed by a unitary time-evolution pUI] , 
recent interest has turned to engineering and controlling 
the time-evolution of open quantum systems. There, the 
goal is to design couplings to an environment, so that 
the many-body time-evolution corresponds to a specified 
Kraus map or, in the Markovian limit, a master equa- 
tion (ME) [5]. In Ref. [5] a scenario of open system dy- 
namics was discussed, where driven-dissipative dynamics 
drives the system of interest into a (desired) entangled 
steady state, and where competition between the coher- 
ent Hamiltonian and the dissipative Liouvillian terms in 
the ME gives rise to a rich nonequilibrium phase dia- 
gram and dynamical phase transitions [7[[S], which do not 
have immediate condensed matter counterparts. Such 
far-from-equilibrium systems may thus open up new per- 
spectives for many-body physics, challenging both theory 
and experiment. In addition, engineered dissipation also 
provides the basis for a dissipative variant of quantum 
computing, as first suggested in Ref. [5]. 

Implementing these concepts of reservoir engineering 
of quantum many-body systems has been so far mainly 
discussed with cold atoms and ions, and first experiments 
have demonstrated driven-dissipative preparation of a 
four-particle Greenberger-Horne-Zeilinger (GHZ) states 
with ions [TU], and Einstein- Podolsky-Roscn (EPR) 
states of atomic ensembles [TTj . In light of the recent 
discussion of cavity arrays and photonic quantum simu- 
lation [T^l |T3] , it seems natural to ask to what extent the 
concepts of quantum reservoir engineering and associated 
nonequilibrium phenomena can be realized with coupled 
photonic systems. In cavity arrays, photon loss appears 
as a natural decoherence mechanism and a steady state 
can only be sustained under nonequilibrium conditions, 
i.e. by continuously driving the system with external laser 
fields O [12] . However, the design of non-trivial cou- 



plings of optical photons to an engineered reservoir, lead- 
ing, for example, to number-conserving or non-local dis- 
sipation processes, is harder to achieve and imposes a 
challenge for these systems. In view of the remarkable 
progress in the field of opto-nanomechanics, we will be 
interested below in setups where the nanomechanical ele- 
ments play the role of quantum reservoirs, and a properly 
designed coupling of nanomechanical oscillators to cavity 
modes results in the desired quantum dissipation. 

In conventional optomechanical systems (OMS) [Ml 
[T7] a single optical mode is coupled via radiation pres- 
sure interactions to the motion of a macroscopic mechan- 
ical resonator, represented, for example, by a moving 
end-mirror or a vibrating membrane |18j inside a Fabry- 
Perot resonator. Over the past years this coupling has 
been successfully employed to cool mechanical systems 
close to the quantum ground state [T9ll23j . using tech- 
niques analogous to laser cooling of atoms. In parallel, 
rapid progress in the fabrication and control of OMS, and 
in particular new designs for micro- and nano- scale de- 
vices [SQlEUGIlGI], have led to a drastic improvement of 
OMS and pave the way for realizing various strongly cou- 
pled [5BH2H] and multi-mode [25H55] scenarios. Here we 
describe the appearance of novel dissipation processes in 
extended OM arrays, where in contrast to OM laser cool- 
ing, now the mechanical systems provide a decoherence 
channel for light. In particular, this mechanism already 
provides a basic building block for reservoir engineering 
for light, and we discuss how this can be used to study the 
dissipation-induced preparation of photonic states with 
long-range coherence |36j . 

The paper is organized as follows. In Sec. [TT] we re- 
view the key concepts behind the engineering of driven- 
dissipative quantum many-body systems. In particular 
we introduce the driven-dissipative Bose-Hubbard dy- 
namics discussed in Ref. [6 for cold atoms, which yields 
a Bosc-Einstcin condensate (BEC) in steady state. In 
Sees. |III| and |IV| we show how the analogous ME for 
driven-dissipative dynamics of photons in an optome- 
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chanical (OM) system can be engineered. In Sec. [V] we 
show that, under continuous incoherent pumping, the 
system undergoes a nonequilibrium phase transition. Fi- 
nally, in Sec. |VI| we summarize our findings. 

II. DRIVEN-DISSIPATIVE QUANTUM 
MANY-BODY SYSTEMS 

In this section we briefly introduce and review the 
many-body ME for driven-dissipative dynamics. The 
specific example is dissipative Bose-Hubbard dynamics, 
as discussed previously for cold atoms in optical lattices 
immersed in an atomic BEC, which plays the role of a 
reservoir of phonons. The dissipative dynamics drives 
the system in steady state into a BEC, representing a 
"dark many-body state." This discussion will set the 
stage for the MEs for photons interacting with an array 
of micromechanical membranes in an optical cavity in the 
following sections. 

A. Engineered steady states 

The time-evolution of a quantum system that is weakly 
coupled to an energy or particle reservoir is governed, in 
the Markov approximation, by the ME 

d tP (t) = -i[H,p(t)}+£[p(t)}, (1) 

where p(t) is the density matrix of the system, % the 
Hamiltonian generating unitary dynamics, and the Li- 
ouvillian £ is a non-unitary, dissipative contribution. 
In general, the Liouvillian C can be written as a sum 
KiA[Kg][p(t)] of terms in Lindblad form 

A[K e ][ P (t)} = 2k t p{t)k\ - k\k tP {t) - p{t)k\k t , (2) 

where {Ki}g is a set of so-called jump operators. In the 
case of a photonic system on a lattice, for example, with 
annihilation operator at on the £th site, local losses of 
photons are described by Kg = hg. The non-unitary part 
of the evolution competes with the interaction Hamilto- 
nian, e.g. the local Kerr nonlinearity H — ¥ £^ a\a\agai, 
and in general the steady state of the system cannot be 
fixed by controlling the eigenstates of % only. Under 
the action of these local jump operators, a pure initial 
many-body state |\E'), extended over several lattice sites, 
evolves into a mixed state. 

Recently, jump operators of a novel kind have been 
designed, which drive the system towards a well-defined, 
pure many-body steady state 

|*)<*| . ( 3 ) 

The engineering of such jump operators is based on 
the combination of coherent driving and dissipation pro- 
cesses, hence the resulting models are properly dubbed 



driven-dissipative systems. The main goal of such reser- 
voir engineering is to obtain many-body states \^>) which 
feature interesting properties, e.g. long-range or topolog- 
ical order. These target states do not necessarily have to 
be associated to the ground state of some physical Hamil- 
tonian and do not result from minimizing energy, but are 
stabilized by a combination of drive and dissipation. 

The case in which |\&) represents an atomic BEC of 
neutral atoms in an optical lattice has been considered 
in Ref. [6 . In the condensed state 

/ l \ N 
|BEC> = ^(^X>!) Ivac), (4) 

all N bosonic atoms occupy the same single-particle 
state, given by the symmetric superposition of the L lat- 
tice sites. In the absence of interactions, the target state 
| BEC) is obtained by the bifocal jump operator 

Kg oc (a\ +1 + a\){a t+ i - a e ) . (5) 

A weak interaction, taken into account perturbatively, 
acts on the steady state as an effective temperature 
fceT oc U and spoils the purity of the final state. When 
the interactions are so strong that they dominate the dy- 
namics, the steady state is mixed and described by a di- 
agonal density matrix. In the intermediate regime where 
the two energy scales compete, it has been found that 
a phase transition takes place pj [5], which shares fea- 
tures of a quantum phase transition in that it is driven 
by the competition of two non-commuting operators, and 
of a statistical phase transition in that the ordered phase 
directly terminates into a strongly mixed state. 

In the case of a photonic system, it is necessary to 
complement the Hamiltonian dynamics with a nonuni- 
tary term C to take into account the decay of the pho- 
tons. More precisely, the lifetime of the photons is usually 
shorter than the timescale oc U~ l over which the unitary 
many-body dynamics takes place. In this case, it is not 
possible to assume that the photons thermalize, or reach 
the ground state of the Hamiltonian T-L, but the complete 
dynamics of the system has to be considered, in the pres- 
ence of pumping and losses. In general, the existence of a 
relation between the ground-state properties of a closed 
many-body system and its open, driven counterpart is 
not obvious, although it may emerge under specifically 
tailored drivings [15l [37] . 

To design a toolbox for the simulation of driven- 
dissipative systems would allow to study the relation 
between the steady-state properties of a nonequilib- 
rium model ([!]) with number-conserving jump operators, 
e.g. Eq. ([5|, and the corresponding open system with 
particle losses. Moreover, the relation between equilib- 
rium and nonequilibrium phases of a model [35] could be 
investigated systematically. More generally, the flexibil- 
ity in the implementation of the jump operators could 
be leveraged to study the effects of the competition be- 
tween different terms in the nonunitary contribution C 
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FIG. 1. Bosonic modes in two bands of a super lattice. Co- 
herent driving of neighboring sites from the ground to the 
upper band and incoherent decay is represented by straight 
and wavy arrows, respectively. The whole array supports a 
chain of A-systems, each defined by a triplet of neighboring 
modes, two of which belong to the ground band. Phase lock- 
ing takes place in the shaded area as a result of the combined 
action of the chain of A-systems. 



to the Liouvillian in a driven-dissipative nonequilibrium 
scenario with stable stationary states, and no immediate 
condensed matter counterpart. 



B. Dark states and A-systems 

The controlled achievement of a given target steady 
state via a certain driven-dissipative scheme is guaran- 
teed in a situation where the stationary state of the dy- 
namics ([I]) coincides with a unique dark state of the Li- 
ouvillian, i.e. for |*) such that K £ \$)=0 for all £ 0139]. 
For example, it was shown in Ref. jB] that the Bose- 
Einstein condensate is a dark state of the jump operator 
can be shown to be unique, and moreover to be 
the unique stationary state of the dissipative dynamics 
(no stationary states other then the dark state exist in 
this case) [39]. The design of jump operators that fea- 
ture a given dark state is thus of central importance to 
the success of a driven-dissipative scheme. Such task has 
been accomplished to generate orbital pairing [3D] as well 
as topological states [41] in the steady state. 

The design of appropriate jump operators is much 
simplified when it can be reduced to the generation of 
quasilocal couplings, which are then repeated periodi- 
cally in the array. For example, the jump operators ([5| 
act pairwise on the lattice sites, annihilating any anti- 
symmetric component in the wavefunction on a pair of 
sites, and mapping it back into the symmetric one. The 
repetition of this quasilocal action results in the phase- 
locking of the wavefunction in the whole array, i.e. the 
generation of the long-range order. 

The jump operator ([5| contains the effective dynam- 
ics of a A-system coupled to a dissipative environment, 
where both the latter and the excited state have been 
adiabatically eliminated. Fig. [T] illustrates the basic pro- 
cess in which bosons in the ground band are coherently 
driven (straight lines) with opposite Rabi frequencies to 
an upper band of excited states. Only the antisymmetric 
superposition between neighboring ground modes is ex- 



cited because of the choice of the Rabi frequencies. The 
excited state decays back to a combination of ground 
states. The spontaneous decay is made possible ener- 
getically by the creation of an excitation in the reser- 
voir, which carries away the released energy (wavy lines) . 
Note that, while energy is exchanged between system 
and reservoir in this way, the system particle number is 
conserved during this process — the system constituents 
merely make transitions between different bands. A de- 
scription of these processes in terms of the jump opera- 
tors ([5]) is valid given two approximations hold: (i) the 
Born-Markov approximation underlying the ME descrip- 
tion requires the separation of the two bands exceeding 
all other energy scales in the problem; (ii) the adiabatic 
elimination of the upper band giving rise to an effective 
single-band description is justified for a driving process 
in a regime where the detuning exceeds the Rabi fre- 
quencies, as well as the additional energy scales involved 
in the dynamics, in particular the many-body couplings 
that arise when several A-systems are joined together. 
We refer to Ref. [5] for a review of the physics of engi- 
neered driven-dissipative systems. 

In Sec. |III| we propose an implementation of the A- 
system using two micromechanical membranes in an op- 



tical cavity, and in Sec. IV we demonstrate that one ob- 
tains a set of jump operators analogous to |5|, once the 
adiabatic elimination has been performed. 



III. OPTOMECHANICAL A-SYSTEM FOR 
PHOTONS 

The implementation of non-trivial dissipative processes 
as, for example, given by Eq. has previously been de- 
scribed in the context of cold atoms by engineering appro- 
priate couplings to a bath of Bogoliubov excitations [6] . 
While in photonic systems dissipation arises more nat- 
urally, it usually appears in the form of single-photon 
decay, and the design of both nonlinear as well as num- 
ber conserving photon processes imposes a challenge. In 
the following, we describe how this can be achieved in 
OM arrays, by making use of nonlinear radiation pres- 
sure interactions between a set of optical and a set of 
damped mechanical modes. 

The basic idea is illustrated in Fig. [2^i for a "mem- 
brane in the middle" setting [18) . where (in a minimal 
version) two semi-transparent membranes separate the 
optical field into two modes 0^2 with frequency u and a 
third mode c with a slightly higher frequency uj' . Oscil- 
lations of the membranes around their equilibrium posi- 
tions then lead to phonon-assisted scattering of photons 
between the modes d\ t 2 and c, whenever the mechanical 
vibration frequency f^M is close to the optical frequency 
difference, J7m ~ — w. When the mechanical systems 
are classically driven this scattering is coherent and — 
in analogy to atomic three level systems — by choosing 
appropriate phases of the driving fields we can select cer- 
tain "bright" and "dark" combinations of a± and 0,2- In 
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addition, the coupling of the mechanical systems to a 
bath of phonons in the support of the membranes, analo- 
gously to the coupling to a low temperature photon bath 
for laser cooling [421 143] , provides a natural dissipation 
channel. Phonon-assisted scattering events followed by 
a phonon decay become irreversible, and thereby intro- 
duce an effective number conserving dissipation process 
for photons. 

In the remainder of this section we describe in more 
detail how the combination of coherent and incoherent 
processes in OM systems can be used to implement a 
driven-dissipative many-body system for photons. For 
this analysis we will focus for concreteness on the "mem- 
brane in the middle" configuration shown in Fig. [2^i and 
consider the extension of this model to a whole OM ar- 



ray in Sec. IV However, our results are not restricted 



to this specific setting and apply to various other OM 
systems. For example, in panels (c) and (d) of Fig. [2] 
we illustrate two equivalent setups, where the OM array 
is implemented using an array of coupled microtoroidal 
cavities or coupled OM crystal cavities (see e.g. Ref. [44] 
and Ref. [21], respectively). Due to the strong confine- 
ment of the optical field on a few micrometer scale, these 
setups enable much stronger OM interactions and can be 
more easily scaled to larger arrays. 



where c and dg are lowering operators for the photonic 
modes localized between the two membranes (with fre- 
quency a/), and between the membranes and the mirrors 
(with frequency cj), respectively. From the very rich spec- 
trum of the system, here we have singled out two frequen- 
cies only, such that the difference 5ui = ui' — cj is almost 
resonant with the membrane frequency CIm, with detun- 
ing Am = 0,ui—8u). We remark here that 5u> is not related 
to the free spectral range of the empty cavity, and the 
equilibrium position of the membranes can be chosen to 
fulfill the resonance condition, exploiting the large num- 
ber of avoided crossings between optical modes [151 143] . 
We assume here that a suitable level structure is found 
in the neighborhood of an avoided crossing, and that the 
rest of the spectrum can be neglected, in analogy to the 
two-level approximation commonly employed in atomic 
physics. 

The frequency of the mode ag depends on the position 
operator Xi, while the frequency of c depends on the 
difference X\ — X 2 . Upon expanding such dependence 
to first order, one finds the linear dispersive coupling 
between the photonic modes and the membranes. The 
resulting Hamiltonian reads 

TioM = -ga\aiX 1 + ga\a 2 X 2 + g'c^c (X 1 - X 2 ) . (9) 



A. The model 

We start out by describing the OM setup that we en- 
vision. We consider two micromechanical membranes 
placed in an optical cavity composed of two macroscopic, 
high-finesse mirrors (see Fig. [2^,). The Hamiltonian of the 
system is of the form 



(6) 



where Hm describes the quantized motion of one vibra- 
tional mode of the membranes, "H p h describes the pho- 
tonic degrees of freedom, and Wom contains the OM 
coupling between the photons and the vibrational mode. 
The position operator Xi for the membranes, with re- 
spect to the equilibrium position, reads = Xyi(b\ + bi) 

in terms of the lowering operators be for the mechanical 
modes, where Xyi = y 'l / '{2m£lyi) . The Hamiltonian 
reads 



M 



(7) 



The membranes can be subjected to periodic driving 
F £ (t) = F(e +int+i w + c.c.) with frequency O, force F, 
and phase (pe- For simplicity, we assume that only the 
phase of the driving is membrane-dependent. The pho- 
tonic part reads 



"H p h = uia\ai + Lua 2 a 2 + o/c^c 



-J{a\c + c'ai + a\c + C* a 2 ) 



(8) 



Finally, the interaction of the vibrational mode of the 
membranes with all the other mechanical modes of the 
structure is modeled by the Liouvillian 

2 

W)] = E + !) A &M<)] + ^ th A[SJ][p(t)]} . 

1=1 

. (10) 

When the mechanical system is in equilibrium with its en- 
vironment, n t h is the occupation of the phononic bath at 
the energy flyiy and 7 is the average damping. However, 
n t h and 7 should be understood as effective parameters 
in the presence of laser cooling and, in particular, 7 can 
be increased and nth reduced below unity [46T - I48] . 

The effect of driving the membranes is to induce tran- 
sitions between photonic states with different frequen- 
cies. To show this, we perform a unitary transformation 
UHU^ of the Hamiltonian, such that the photonic fields, 
to linear order in J/6u>, transform as 

CL\ 1 — ^ Q/\ — (J/(5w)c, a 2 1— > a 2 — {J/8ui)c, 

c4 — c— (J/5u])ai- (J/5w)a 2 . (11) 

(In the following, only the transformed operators will be 
used so, to avoid cumbersome notation, we do not intro- 
duce new symbols.) The transformation has been chosen 
to diagonalize the hopping part of the Hamiltonian ^ 
that, in terms of the new operators, reads 

"Hph = Loalax + ua\a 2 + u'tfc + 0[(J/Slo) 2 } . (12) 

The quadratic part in J/Su>, omitted from the Hamil- 
tonian for the moment, contains hopping between the 
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ilar approach has been considered in Ref. [3D]. A more 
rigorous derivation, starting from an ad hoc quantization 
of the photonic modes in the presence of a moving mem- 
brane, has been provided in Ref. |51j . 

It is now convenient to perform a displacement be 1— ¥ 



-iflt 



be of the lowering operators of the me- 
chanical mode, where Xe = e _l<w -FX|f/(n — Qjvr+ry/2). 
The displacements removes the driving term from the 
Hamiltonian ([7]) , in the rotating- wave approximation. In 
terms of the displaced fields, the Liouvillian ( 10 ) is in- 
variant in form, while the OM Hamiltonian ( 13 1 can be 
decomposed as Hou = T~Lt> + %vib, where 



£e' int c\a 2 - ax) + H.c. 



(14) 



with the effective driving strength £ = —gJXi/du, and 



{J/8uj)X u [c\g'a\ + (g + g')a\]b\ 
-c[g'a{ + (g + g')ai}bl+Ji.c.} . 



(15) 



FIG. 2. a) OM setup for implementing engineered dissipa- 
tion processes for photons. Three optical modes ai,2 and c are 
delimited by two membranes (vertical thick lines) at positions 
Xl and X2 and the two fixed end mirrors of the optical cavity. 
The modes are tunnel-coupled with amplitude J and their fre- 
quencies are modulated by oscillations of the membranes with 
frequency Qm around their equilibrium positions, b) The re- 
sulting energy level diagram assuming degenerate modes a\ 
and a,2 and a slightly higher frequency uJ ~ u + Q,m for the 
intermediate mode c (see text for more details), c) and d) Al- 
ternative implementations using an array of coupled microto- 
riodal cavities or coupled OM crystal cavities, respectively. In 
both cases the optical fields are confined inside the structures, 
coupled among each other via their evanescent fields and in- 
teract with localized mechanical vibrations of the structures 
via optical gradient forces. 



degenerate modes a#. The OM coupling ([£]) reads 

Hom = -ga\aiXi + ga\a 2 X 2 + g'$c (Xl - X 2 ) 

+(J/6uj)c i [g'a 2 X 1 - g' a x X 2 + (g + g')«iXi 
-{g + g')a 2 X 2 ] + H.c. + 0[{J/5lo) 2 ] . (13) 



In the following, we neglect the first line of ([13]), which 
contains non-resonant terms that account only for a 
renormalization of the bare frequencies in Eq. (12 1 (see 



also Ref. |49j). The remaining lines, on the contrary, rep- 
resent a "three- wave mixing" effect, in which the energy 
of a phonon is absorbed (emitted) to scatter a photon to 
a state of higher (lower) energy. In the present deriva- 
tion, the phonon-mediated coupling between modes with 
different energy arises because of the mode mixing due 
to the hopping terms in the original Hamiltonian. A sim- 



The first part describes photonic transitions induced by 
the membrane driving, while the second part contains the 
coupling between the photons and the quantized vibra- 
tional mode. In Eq. ( 14 1 we also made the specific choice 



<Pi = 0, which means that the forced vibration is actuated 
with the same phase on both membranes. We see that, 



with this choice of phases, a state |^) cx (a\ + a 2 ) |vac) 
[cfr. Eq. B], where all the photonic population is con- 
tained in the mode given by the symmetric combination 
of the modes ai and a 2 , is unaffected by the Hamilto- 
nian %d- In this sense, the symmetric state is dark 
with respect to the mechanical driving, which excites 
only photons in the antisymmetric mode to the upper 
mode c. Finally, it is convenient to remove the rotating 
phase from the Hamiltonian ( 14 ) using a new frame for 



the p hoto nic operators, such that the photonic Hamilto- 
nian (12 1 reads "Hph = — Ac^c, with A = CI — 5uj, and 
V. D = £$(a 2 — ax) + H.c. 

Resonant "three-wave-mixing" processes in OM sys- 
tems have been previously considered for the transduc- 
tion of photons and phonons [52] and have been pro- 
posed for enhancing nonlinear effects in strongly cou- 
pled OM systems [32l [33]. Here we are interested in 
the regime where the mechanical modes, which are in- 
volved in this process, are strongly damped and there- 
fore act as a reservoir for the photons. Assuming that 
the coupling strength between the photonic and phononic 
modes is smaller than the decay width of the phonons, 
i.e. J\g\/8ui, J\g'\/8uo <C 7, we perform the Born-Markov 
approximation and trace out the phononic modes be. As 
a result the reduced density matrix p p h{t) = Try[[p(t)] 
for the photons obeys the ME |l]) with the Hamiltonian 



G 



H = H p h + tiu and the Liouvillian C = £ c ff, where 

2 

C cB [p P h(t)} = £)f (nth + l)A[i^][p ph (*)] 

2 

+ ^fn th A[£j][p ph (t)] , 



(16) 



with the jump operators 



K x =c[al + {l+g/g')a\] 
K 2 = c[a\ + (l + g/g')a 2 } 



(17) 



The effective damping rate in the Liouvillian ( 16 ) reads 
r ~ 2(g' XmJ/Suj) 2 /j. We point out that our derivation 
holds at leading order in J/8u>, i.e. linear for the Hamil- 
tonian and quadratic for the effective Liouvillian (16 1. 



Higher-order corrections to the structure of the coupling 
Hamiltonian would yield corrections to orders higher 
than two to the Liouvillian, and hence can be neglected. 
We also note that, in the limit g' 3> g, the two jump 
operators are equal and, after the adiabatic elimination 
of the mode c (see Sec. IIIBl, they reproduce the form 
^ that was studied in the atomic context [BJ. 

Finally, coherent excitations of the cavity modes by 
external laser fields can be described by 



n co = y^£i(t)q| +H.c. , 



(18) 



with ee(t) = eie~ lQjlt in the rotating- wave approxima- 
tion. Losses of photons through the mirrors of the cavity 
are described by the Liouvillian 

2 

£ioss[p P hW] = £«A[a/][pp h (t)] + KA[c}[ Pph (t)} , (19) 
e=i 

where we assume, for simplicity, that the decay width n 
is the same for all modes. 

In the following, we omit the subscript in the photonic 
density matrix and redefine p p h — > p, which all together 
obeys the ME 

p = -i[tiph +Hd+ tico, P] + £e$P + AossP- (20) 

It describes the open system dynamics of the two cavity 
modes Si and S2, which are externally driven and coupled 
dissipatively via excitation and successive decay of the 
intermediate mode c. 



B. Adiabatic elimination of the excited mode 



The full model for the optical modes given in Eq. ( 20 ) 



can be further simplified if the population in the excited 
mode c is negligibly small. This is indeed the case for 
appropriate choices of parameters, as we point out in the 
following sections. In this case, we may proceed to the 



adiabatic elimination of the mode c to obtain an effective 
ME for the modes di and S2 only. 

We start out by introducing the projection super- 
operators [S3] Vp = Tr c [p] <8) |0) c (0| onto the vacuum 
state of mode c and its complement Q = 1 — V. By ne- 
glecting the weak external driving and the cavity losses 
(which just lead to an overall decay of all modes) in the 
present derivation, the equations of motion (EOM) for 
the projections read 

d t Vp = VC B Qp + VC cS Qp , 

d t Qp = Q(£ ph + C cS + Cb)Qp + QC B Vp , (21) 

where Cd and £ p h are the Liouvillians corresponding to 
Hu and H p h, respectively. For simplicity we consider 
the limit g' g and n t h = 0. Note that here the V 
and Q subspaces are coupled both by the coherent linear 
driving in £d and the nonlinear decay in C e s- Our goal 
is to solve the EOM for Qp perturbatively in the driving 
strength £, without assuming that the decay is small. 
To this end, we decompose Qp(t) = Qi(t) + Qvif) in a 
first- and second-order contribution, which are obtained 



by formally integrating Eq. (21 ) and read 



Qi(t)= f dTe {c ^ +c ° !i)(t - T) QCvV P {t), 

Q 2 (t) - [ dre^+^^-^CuQ^T) . (22) 
J — 00 

These expressions are inserted into the equation of mo- 
tion for V p and the exponentials are explicitly evaluated 
using the decomposition 



(23) 



where the indices s and a refer to the symmetric a s = 
(fii + a2)/v / 2 and anti-symmetric a a = (Si — S2) /v2 
mode, respectively. As a result we obtain an effective 
ME for the reduced density operator of the S modes, 
p a = Tr{Vp}, which can be written in the form 



2£ 2 A 
\ 2 +f 2 



ala a ,p a 



+ 8£ 2 fA[K}p a + £ lossPa , (24) 



where f = 4f (S|S S + 1) and we have introduced the jump 
operator 



K 



1 



iA + 4r(SjS s + l) 



(25) 



In the limit A 3> Y(a\a s ), the Hamiltonian reduces to 
a shift of the antisymmetric mode, while the jump op- 
erator recovers the form i.e. K cx SjS a . For larger 
r, additional higher order nonlinearities arise, which, for 
example, can lead to pumping-induced anti-bunching ef- 
fects, as described in Ref. [36] . 

In currently realized OM devices with membranes in- 
side a Fabry-Perot cavity [TS] [3S] the optical free spectral 
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range is in the GHz range and flexural mechanical modes 
have frequencies Qm — MHz. At the same time, very low 
cavity decay rates n ~ 50 kHz <C can be achieved, 
such that the separation of frequency scales assumed in 
our derivation can be realized. For typical mechanical 
quality factors of Q ~ 10 6 , the bare mechanical decay 
rate 70 can be as low as 1 Hz, but it can be tuned up to 
7 ~ f^M by applying additional laser cooling SSI 27] or 
feedback damping [381 E3] • This would also reduce the 
effective temperature to mean occupations nth — 0(1). 
Note that although in these systems the OM coupling 
is rather low, g ~ 1 — 10 Hz and the resulting dissi- 
pation r ~ g 2 /l "C n, the nonlinear effects we are in- 
terested in will depend on the ratio r(aJa s )/K, which 
can exceed unity for larger photon numbers. Further- 
more, in several solid state realizations of OM devices, 
such as microtoroidal cavities or OM crystals [see Fig. [2] 
(c) and (d)], OM coupling constants g of few kHz up to 
~ 1 MHz [2UJ [3TJ [25] have been demonstrated, and even 
the single-photon strong-coupling regime g > k is within 
reach of optimized devices. In such devices, the regime 
r ~ k is in principle accessible, while at the same time 
the high mechanical frequencies SIm — GHz [3TJ [25] lead 
to equilibrium occupation numbers n t h < 1 at cryogenic 
temperatures. This parameter regime will be considered 
in most of our numerical simulations below. 



C. Numerical results 

We study now in detail the building block of the 
many-body driven-dissipative system, i.e. the effective 
A system for photons introduced in Sec. |III A| accord- 
ing to Eq. (20 1. First we consider the case in which 
the lower mode Si is coherently driven. The coher- 



ent pumping is resonant with the single mode, which 
means that it is equally detuned from both the symmet- 
ric (Si + 02) /v^ and the antisymmetric (Si — a%) /\/2 
modes, which are split by 2J. The splitting arises 
from the residual coupling in Eq. (12), which we write 

"Hhopp = -J(a\a 2 + H.c). 

The EOM are solved using a standard fourth-order 
Runge-Kutta method, truncating the Fock space of each 
mode to 4-6 states. The coherent pumping couples 
equally to the symmetric and antisymmetric mode, how- 
ever, because of the OM driving and dissipation, the 
steady-state population in the symmetric mode is sub- 
stantially larger than the population in the antisymmet- 
ric mode. During the time-evolution, the population in 
the upper mode c can be made negligibly small, employ- 
ing a large detuning A of the OM driving. 

In Fig. [3] we present a clear signature of the popula- 
tion imbalance obtained by means of the OM driving. 
More precisely, we compute first the connected two-time 
Green's function 

g^\tM) = ([4w - <4(*))p2(*o) - <s 2 (i ))]) (26) 

in the mode S2, where the average is computed on the 




FIG. 3. Spectrum of the light radiated from the &2 mode. 
The ai mode is coherently pumped with strength ej = 0.1T, 
all the modes decay with amplitude k, — 0.1F. The OM driv- 
ing is £ = 0.5f (solid line) O.Olf (dashed line), with J = 0.2F. 
The system is first evolved to the steady state at t^F = 50.0 
and subsequently to tT = 500.0 to compute the two-times 
correlation functions with sufficient spectral resolution. 



steady state at t = to. The absolute value of the Fourier 
transform with respect to the time t is then normalized 
to unity to give the cavity output spectrum S(uS). We 
see that, in the absence of OM driving, the spectrum is 
composed of two equal peaks which represent the sym- 
metric and antisymmetric mode (at lower and higher fre- 
quency, respectively). Switching the OM coupling on, 
the antisymmetric peak is strongly suppressed, due to 
the efficient scattering of photons into the symmetric 
mode. The suppression of the peak corresponding to the 
antisymmetric mode is an unambiguous consequence of 
the phase-locking between the two sites and clearly dis- 
tinguishes our setup from a more conventional situation 
where two independent cavities are pumped coherently, 
without being phase-locked to each other. The effect of 
switching on the OM coupling is studied in more detail in 
App. [Bj where the entanglement between the symmetric 
mode and the rest of the system is discussed as well. 

While in Fig. [3] we show the noise spectrum on one 
side of the system while both symmetric and antisym- 
metric modes are pumped coherently, in Fig. [4] we il- 
lustrate the effect of pumping the antisymmetric mode 
only, with variable detuning. We use opposite Rabi fre- 
quencies in the two modes 1=1 and I = 2, so that 
the coupling to the symmetric mode vanishes. The pop- 
ulation in the steady state shows broad resonances cen- 
tered around the frequency J of the antisymmetric mode. 
While the broadening of the noise spectrum in Fig. [3j at 
zero temperature, is determined by the photonic decay 
in each mode, the broadening in Fig. [4] is given by the 
larger scale T. It is interesting to note that this proce- 
dure shows a resonance in the symmetric mode which is 
centered around the energy of the antisymmetric mode 
as well. This is due to the fact that the excess energy 
is taken away by the phonons that are emitted following 
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FIG. 4. Main panel: steady-state populations ng for the sym- 
metric (solid line), antisymmetric (short dashed), and excited 
(long dashed) mode, in the presence of coherent pumping of 
the antisymmetric mode with detuning Ai. Inset: peak value 
of ni at detuning Ai — J — 2.0T with finite thermal oc- 
cupation n t h of the phononic reservoir. The strength of the 
coherent pumping is eg = 2.0k, with £ — 2. OF and F = 10.0k. 



vious case, and only the final results are reported. The 
time-evolution of the extended system, however, is much 
more complex, because the photonic field is inhomoge- 
neous. The dissipative mechanism must single out the 
symmetric mode in the lattice [with zero quasimomen- 
tum k — 0, cfr. Eq. Q] among the many modes of the 
Brillouin zone, while only the central (symmetric) and 
edge (antisymmetric) Brillouin-zone modes are present 
in the two-membrane system. We show numerically that 
the driven-dissipative dynamics achieves phase-locking of 
the photonic field throughout the array, i.e. the photonic 
field converges towards the k = state. 

The membranes are disposed in a regular array of pairs 
at positions X 2 t-\ and X 2 g (see Fig. [5]) and define two 
sets of localized photonic modes (ag and eg) with differ- 
ent frequencies (u and u' , respectively). Because of pho- 
ton transmission through the membranes, the spectrum 
splits into two bands centered at the two frequencies u> 
and ui' . The splitting depends on the hopping amplitude 
J between the modes, and we assume that it can be ne- 
glected with respect to the difference Sui = ui' — u>. The 
transformation (111 now reads in general 



/ 
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X 2 e-i X 2 g 



FIG. 5. Schematic representation of the system under consid- 
eration. An array of membranes (thick vertical lines) , grouped 
into doublets at positions X21-1 and X21, within an optical 
cavity, define two sets of quantized localized light modes ag 
and eg. In Sec. |III| we discuss the case in which only two 
membranes are present. 



the (virtual) excitation of the mode c. The latter mode, 
as is apparent from the figure, is always negligibly popu- 
lated. Finally, the effect of the finite temperature of the 
phononic reservoir is investigated in the inset of Fig. [4] 
We see that the driven-dissipative mechanism is robust 
in the presence of low thermal occupations n t h, as the 
magnitude of the peak corresponding to the symmetric 
state remains substantial. 



IV. MEMBRANES ARRAY 

In this section we discuss the extension of the A- 
system for photons of Sec. |III| to the case of an array 
of 2L membranes. The model for the membranes ar- 
ray [34 , 35j is a straightforward generalization of the pre- 



ag H> wen + rag + wdg, eg i-> w ag + r eg + w ag + \ , 

(27) 

and the following results hold to the first nonzero order 
in w and w' . A specific example of the transformation is 
given in App. |Xj The photonic part of the Hamiltonian, 
in the rotating frame, reads ?i p h = — ^2g=± AcJq, while 
the driving part (flil) becomes 



L-l 

E 

1=1 



[{-Y£c\{c 



ag+i) 



H.c. 



(28) 



with £ — 2gXw*r. Here we assume that each pair of 
membranes vibrates rigidly together, while neighboring 
pairs have opposite velocities. More precisely, the phase 
of the driving is f-u-x — f2i = if £ is odd, and ip 2 g-i = 
^21 — if t is even. In the effective Liouvillian ( 16 ), the 



summations extend to I = L with the jump operators 

K 2 g-i = xa\(ct-i + eg) + (a\ + a) +l )cg, 
K 2 g = xa\ +1 {cg + cg + i) + (a] + a\ +l )cg , 

(29) 

with x = gw*r I '(g'w'r'*) and the effective damping T ~ 
2|s'AWrf/ 7 . 

We first investigate the time-evolution of an inhomo- 
geneous initial state, in the absence of pumping and dis- 
sipation, in the semiclassical approximation. To derive 
the semiclassical EOM, we first compute the EOM for 
the averages ag = (ag) and (eg) and we substitute each 
field operator with the corresponding c-number. The so- 
lution of the EOM in Fig. [6] shows that the system is 
driven towards a state where the photonic occupation in 
the ag modes is homogeneous, while the occupation in 
the eg modes decreases exponentially. The relative phase 
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2 3 4 5 
ff xlO 4 



2 3 4 5 
ff xlO 4 



FIG. 6. Time-evolution of the photonic population in the ex- 
cited [panel (a)] and lower [panel (b)] levels, and of the relative 
phase [panel (c)] between the semiclassical field on neighbor- 
ing ground sites, computed according to the average-number- 
conserving semiclassical EOM. L = 12 sites (corresponding to 
2L membranes) are used with periodic boundary conditions. 
The OM driving is 0.2f with detuning -10. Of . 



<pi = Arg[o^o^ +1 ] between neighboring ai modes van- 
ishes, demonstrating the onset of phase-locking in the 
array and the generation of the state with quasimomen- 
tum k = 0. 

Given the exponentially suppressed occupation of the 
modes q, visible from the solution of the semiclassical 
EOM and from the exact diagonalization discussed in 
Sec. |III[ we may proceed to the adiabatic elimination 
of such modes. This allows us in Sec. [V] to focus on the 
many-body quantum dynamics of the modes ag only. The 
EOM for the operator eg close to the the steady state is 

d t h(t) ~ *Act(t) - i(-YS[a e (t) - a e+1 (t)} . (30) 

Here we neglect the contribution of the dissipation T with 
respect to the detuning A, as demonstrated in Sec. |IIIB| 
in the case of two membranes. We may hence substitute 
the steady-state form 



(31) 



into the jump operators (J29J). The adiabatic approxima- 
tion correctly captures that the average (ce(t)) vanishes 
when on — ctt+i, i- e - on the steady state, as shown by the 
solution of the complete semiclassical equations in Fig. [6] 
The effective jump operators for the modes ag read 



K- 



21-1 



£ 



K 



2C 



H l - 



aldg 



(2x + l)at , ,o w 



<X 



At 



At 



a e ae + i + x a g+\ a i+2 



(32) 



We notice that these jump operators are considerably 
more involved than ^ as they couple up to three neigh- 
boring sites each. Each jump operator Kg describes the 
dissipation arising from the phonons in the £th mem- 
brane. Dissipation processes due to different membranes 
are of course independent, and are then described by 
different Lindblad terms The existence of a dark 

state when nth = (i- e - when the temperature of the 
phononic reservoir vanishes) can be easily seen in mo- 
mentum space, where the jump operators read (apart 
from an irrelevant overall sign) 

ku-i = ^4 afe( _ 1) ^ e ^-D e -^+i) 



x e 



K- 



q.k 

*(*!+?) n 

q,k 

x[l + X 



x) 



Ak 



e iq xW k 



2 , = ^E« 



e* k X + e^xW k -l} , 



(33) 



with ag = ^2 e +lqe a q /y/L. Indeed, we see that the state 
with momentum k = 0, which is the analogous in the 
lattice of the symmetrized mode (Si + a?)/ y/2 considered 
in Sec. |III| belongs to the kernel of the jump operator, 
and hence its population is never depleted by the driven- 
dissipative dynamics. 



V. DYNAMICAL PHASE TRANSITION 

In this section we consider the dynamics of the mem- 
brane array, introduced above, under continuous inco- 
herent pumping and photonic losses. In the steady state, 
the light intensity is governed by the rate of pumping 
and losses only. However, the light coherence between 
different sites depends dramatically on the strength of 
the OM driving and vanishes below a critical point. We 
argue that this effect can indeed be interpreted as a dy- 
namical phase transition, which can be used to probe 
the effectiveness of the driven-dissipative mechanism in 
a strongly nonequilibrium setup. 

We first discuss the Gutzwiller approximation, which is 
used to tackle the highly nontrivial problem of the time- 
evolution of the many-body system. Then, we present 
the numerical solution to the EOM and, finally, a com- 
pact analytical method is presented, which yields a global 
picture for the the steady-state light coherence proper- 
ties. 



A. Mean-field approximation 



To investigate the quantum dynamics induced by the 



jump operators (32) beyond the semiclassical approxima- 
tion used in Sec. IV we resort to a mean-field approxima- 
tion where the density matrix p is factorized in real space 
into a product p ~ pi of reduced local density matri- 
ces pi = Tr^eP on the £th site [3[S]. We remark that this 
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approximation goes beyond the semiclassical approxima- 
tion, as it conserves the full quantum structure of the lo- 
cal Fock space, so that local correlation functions are not 
factorized. This method has been successfully employed 
in Refs. [3 [S] to study the time-evolution of number- 
conserving driven-dissipative systems. The hypotheses 
that support the usage of the mean-field decoupling are 
valid in the present case as well, because the addition 
of local pumping and dissipation to the theory does not 
crucially rely on the representation of the nonlocal cor- 
relation functions, which are severely approximated in 
this mean-field scheme. The equation of motion for the 
reduced density matrix pt(t) follows from 



dtpi{t) = Tr#[£[p(t)]] = C t \p e (t)] 



(34) 



using the expression detailed in App. [D] In the following 
notation, the index of the site can be omitted from the 
operators, because the problem has been reduced to a col- 
lection of (nonlinear) local problems. We also presently 
reduce to the homogeneous case, i.e. we assume that the 
average {■ ■ - )i on the tth site is independent of I. The 
mean-field approximation applied to the Liouvillian with 
jump operators (32) gives 



[r(n th + l)i 



(i) 

ij 



Tn th L 



(2) 



-P(*)}] 



(35) 



with T = T£ 2 /A 2 and A = (1, a, a\ a) a). In the fol- 
lowing, we focus on the zero-temperature limit nth = 0. 
The 4x4 matrices L^ 1,2 ) are reported in App. [dJ They 
depend on the expectation value of several (at most cu- 
bic) products of creation and annihilation operators and 
are functions of the density matrix p. The equation of 
motion for p is hence nonlinear in the mean-field approx- 
imation. For n t h = 0, the existence of the dark state 
with vanishing quasimomentum k = is reproduced due 
to the fact that the Liouvillian vanishes when applied to 
the pure coherent state. 

The photonic modes a are subjected to local losses 
at rate k, which reduce the total particle number. To 
prevent the photonic population to vanish in the steady 
state, photons are injected into the system at rate T by 
means of an incoherent pumping scheme. The Liouvillian 
operators for such processes read 



pump 

[p(t)] = TA[dJ][p(t)}, 



AossIpW] = «A[o][p(t)] 



(36) 

In the following, we consider the full dynamics in the ar- 
ray in the presence of the number-conserving Liouvillian 
(35) and of the pump and loss terms (36). The latter 



terms induce a continuous flow of photons through the 
system, which has a substantial impact on the properties 
of the steady state. Indeed, the effect of the Liouvil- 
lian (35) goes towards the creation of a state featuring 



long-range phase locking in the array, while the pump 
and loss terms induce local dephasing. The system ex- 
hibits then a competition between the two contributions 



j i i i_ 




5 10 15 20 25 30 5 10 15 20 25 30 

t K t K 

FIG. 7. Time-evolution of the coherent fraction \a\ 2 jn (left 
panel) and second order correlation function </ 2 ' at zero delay 
(right panel), with incoherent pumping T = 0.5k, for three 
different values of the rescaled driving strength F' ~ 0.5k 
(long dashed), F' ~ 1.5ft (short dashed), and F' = 2.5ft (solid 
line), at zero temperature nth = 0.0. The system is perturbed 
for tn < 1.0 by a weak coherent pumping field. 



to the non- unitary dynamics 9J. Such competition re- 
sults in a dynamical phase transition, as it is shown in 
the next section. 



B. Onset of coherence in the presence of 
incoherent pumping 

In the absence of the dissipative OM coupling T = 0, 
the combined effect of incoherent pumping and photon 
losses is to create a thermal state [55] with steady-state 
density = T/(k — T), where we require T < k for 
the system to be stable. The density matrix of the ther- 
mal state is diagonal and the average a = (a) vanishes, 
indicating the absence of coherence. The "coherent frac- 
tion" \a\ 2 /n achieves the maximal value 1 in the coherent 
state, which is assumed in the semiclassical approxima- 
tion. In the Gutzwiller mean-field approximation for the 
density matrix, both these limiting cases are contained: 
the Fock space onsite density matrix can describe both a 
thermal state with matrix elements p n , n ' = Pn^nn'i Pn = 
n n /(n + l) n+1 as well as a coherent state p = |4 r ) C oh(^ , |, 

l*)coh = Er=o e ~' a|2/2 («™/^)l^) indicating the pres- 
ence of off-diagonal order. The coherent fraction mea- 
sures the long-range coherence in the array, because the 
correlation functions factor as (al,ai) ~ (a^)* (ag) on 
the coherent state. In this sense, the ratio |a| 2 /ro corre- 
sponds to the condensate fraction in the mean-field the- 
ory of weakly interacting Bose-Einstein condensates (see 
e.g. Ref. [55]). 

The main result of the present section is that a crit- 
ical value T' crit ~ k exists, above which the condensate 
fraction is finite. (The rescaled driving strength r is in- 
troduced in App.[C|) The phase-locking effect engineered 
into the Liouvillian ( 35 ) overcomes the phase randomiza- 
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FIG. 8. Dependence of the coherent fraction \a\ 2 /n on the 
rescaled driving strength V' , at zero temperature nth = 0.0, 
for T = 0.5k. The system is initialized in the vacuum state, 
perturbed for te < 1.0 by a weak coherent pumping field. 
Stroboscopic plots of the profile at increasing times (dots, 
dash-dots, short dashes, long dashes, solid line) are shown to 
converge towards a steady-state profile. 



tion induced by the photon losses and the pumping ( 36 1 , 



generating a state with long-range order in the array. 
This result is first illustrated in Fig.[7J where we show the 
time-evolution of the coherent fraction for several values 
of the coupling V . We also show the second-order corre- 
lation function at zero delay g^ 2 \r — 0) = {a*a)aa)/n 2 
for the same parameters. As the coherent fraction in- 
creases, g( 2 ' decreases from 2 towards 1, signaling a tran- 
sition from a thermal to a coherent state. The steady- 
state value of the coherent fraction as a function of the 
driving strength T' is finally shown in Fig. [8j which can 
be seen as a nonequilibrium phase diagram in one vari- 
able. The incoherent and coherent phases are defined by 
the "order parameter" a being zero and nonzero, respec- 
tively. 

We remark that the transition arises from phase- 
locking, which is a quasilocal, single-particle phe- 
nomenon, and not from coupling some local degree of 
freedom to a collective (mean-field) excitation. This con- 
sideration sets the present setup apart from other exam- 
ples [571 HH] of spontaneous onset of an order parameter 
in a many-body system. 



C. Analytical interpolation with displaced thermal 
states 

To provide an analytical description of the nonequi- 
librium phase transition presented above, we derive now 
an effective dynamics for the field expectation value a. 
The effective dynamics that we construct is such that a 
converges to a — in the normal phase, and to \a\ > 
in the coherent phase. Focusing on the parameter region 
around the critical point, a performs an overdamped mo- 
tion on a energy potential U(a) ~ C2|a| 2 + C4|a| 4 . Our 



totic dynamics of the system in terms of its collective 
variables. In this sense, a plays the role of the order pa- 
rameter of the dynamical phase transition. In analogy 
to Landau theory, where the minimization of the free en- 
ergy gives the equilibrium phase, the fixed point of the 
effective dynamics gives the steady state of the nonequi- 
librium system. 

To construct the effective dynamics, we make the 
ansatz that, close to the steady state, the density ma- 
trix p depends only on a and can be approximated by 
the form 



p(a) ee V{a)p{n)'b{-a) 



(37) 



where the action of the Glauber operator T> is 
T>{— a)aD{a) — a + a, and p is the thermal density ma- 
trix with average density n. The total photonic density 
n — n + \a\ 2 features an incoherent and a coherent con- 
tribution (that correspond to the thermal and condensed 
fractions, if the same formalism is adapted to describe ul- 
tracold bosonic matter). Particular cases of this ansatz 
include the coherent state, where n = |a| 2 , and the ther- 
mal state where a = and the density matrix is diagonal. 

The time-evolution dtn(t) = — 2(n — Y)n(f) + 2T of the 
total density is determined by the incoherent rates of in- 
jection and leakage of photons and gives = T/(k — T) 
in the steady state. In this respect, the incoherent part 
of the time-evolution plays a role analogous to the chem- 
ical potential in a system of massive bosons at equilib- 
rium. The fact that the total number of particles does 
not depend on the driving strength T suggests that we 
may define the effective dynamics (evolving in time r) in 
such a way that the thermal state is always the instan- 
taneous solution of the EOM and d T n = (see details 
in App.[C]). In the proximity of the critical point, where 
the time-evolution of the order parameter a is expected 
to be much slower than any other macroscopic degree of 
freedom of the system, this is certainly a good approx- 
imation. In the remaining part of the phase diagram, 
this approximation can be seen as a coarse graining of 
the dynamics. Using the value n = T/(n — T + |a| 2 r'), 
which fulfills the ansatz in the coarse-grained approxi- 
mation (see App. [C|), the effective equation of motion for 
the order parameter reads 



Or 



-2(«-T) 



T 



T 



T + 



2p 



(38) 

The stability of the system (meaning that the coherent 
density does not increase arbitrarily) is granted by k > T 
because i) the prefactor multiplying the square brackets 
is negative and ii) the second term in the square brackets 
vanishes as \a\ 2 — > oo, with the sum of the two remaining 
terms positive. It is interesting to see that, while the 
stability of the system for r' = is obvious, for T' > 
the effect is nonperturbative, as it stems from a term 
cx \a\~ 2 and cannot be obtained from a generic expansion 
in the Landau form. The order parameter around the 



construction offers an effective description of the asymp- critical point follows the overdamped equation — d 2 \a\ 
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-9 T |a| — dU /d\a\, for m — > oo, with the effective energy 



U =(k-T) 1 -^- + T j 



|a|_ 
4 



(39) 



The solution for the steady states of the effective dy- 
namics (which coincide by construction with the min- 
ima of the effective potential U) then reads \a\ 2 = for 
V < T' crit and H 2 = noo(l - T' CT JT') for V > T' cl , t , 
with T' clit = (k — T) 2 /T. It is interesting to note that, 
for any k, one can drive the system strong enough to 
overcome thermalization and enter the phase where the 



transition takes place. Eq. (38) is applicable to all values 



of the driving V , features the correct results in the limits 
r — > and r — > oo, and correctly captures the existence 
of a nonequilibrium phase transition in the model. The 
profile of | a | 2 obtained by the numerical solution in the 
close vicinity 6V / k < 0.1 of r£. rit , though, is more com- 
plex than the linear slope which follows from Eq. ( 39 ) . 



VI. CONCLUSIONS AND OUTLOOK 

In conclusion, we have devised a scheme to implement a 
driven-dissipative dynamics for photons in an OM setup. 
The mechanical elements offer both a way to drive the 
photons between cavity modes, and a natural reservoir 
where energy can be dissipated. Due to the huge energy 
difference of the photonic and the phononic modes, the 
engineered effective dynamics for the photons is number- 
conserving. The basic block of our scheme is an effective 
A-system for photons, which generates a dark state con- 
sisting of the symmetric superposition of two photonic 
modes. Extension of the scheme to an array of mem- 
branes is shown to generate a state with long-range coher- 
ence, where the phase-locking is achieved throughout the 
array. The interplay of the engineered driven-dissipative 
dynamics with the losses of photons out of the cavity 
and incoherent driving of the photonic modes results in 
a dynamical phase transition. We provide a compact 
analytical description of such transition, in terms of an 
effective potential that is minimized in the steady state. 

We have focused on an implementation based on mi- 
cromechanical membranes in an optical cavity, where the 
photonic modes are delimited by the membranes and the 
mirrors of a Fabry-Perot resonator. This kind of setup 
offers extremely low optical losses and large tunability in 
the position of the mechanical oscillators. It is an ideal 
candidate to test the building block of our setup, which 
consists of two membranes in a cavity. However, other 
OM implementations, such as arrays of nanotoroids or 
OM crystals, provide a better scalability, and might be 
more suited to investigate the onset of long-range coher- 
ence in the array. The scheme that we devised is generic 
and could be extended to such implementations as well. 

OM setups are complementary to other quantum simu- 
lation frameworks, such as cold atoms or circuit-QED, as 
to which kind of many-body phenomena can be realized. 
The possibility to manipulate the energy spectrum by 



geometrically displacing the mechanical oscillators, and 
the naturally arising coupling to well-controlled phononic 
reservoirs, make OM setups ideally suited to implement 
nonequilibrium models. In this work we considered 
the time-evolution and stationary state of a many-body 
model, in the presence of continuous flow of particles 
through the system. The existence of a dynamical phase 
transition, which is not present in the number-conserving 
dynamics, suggests that it is interesting to study other 
many-body models under such conditions as well. Fi- 
nally, once strong-coupling of photons and phonons is 
available in OM devices, the effect of dissipation on well- 
known equilibrium quantum phase transitions can also 
be investigated along these lines. 
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Appendix A: Localized modes in the array 

In this appendix we substantiate the existence of the 
two families of modes hi and q discussed in Sec. |IV[ The 
procedure that we use relies on the possibility to define 
decaying modes a for light between leaking mirrors, as 
detailed in Ref. [SS] . The existence of hopping processes 
between the internal modes of a cavity and all the ex- 
ternal modes of the electromagnetic spectrum, even at 
different frequencies, is commonly used to describe leak- 
ing cavities from the quantum noise perspective [551 160) . 
In our case, with the hopping term we model the exis- 
tence of nonzero matrix element between decaying modes 
bounded by membranes and mirrors, with different en- 
ergy. In the absence of external driving, these transitions 
are of course only virtual and the real transport through 
the array takes place via double-tunneling events of am- 
plitude J 2 . 

To discuss the form of the modes, we use a single- 
particle picture with Hamiltonian 



H = 



I 


-1/2 


J 









J 


+1/2 


J 


n 


\ 





J 


-1/2 




J 





J 


+1/2/ 



(Al) 



which describes an array of double wells, where the en- 
ergy difference Suj between the wells has been taken equal 
to 1, and J is the hopping between the wells in units of 
Slu. The basis in which the Hamiltonian is written sup- 
ports the definition of the modes ag (with energy —1/2) 
and C£ (with energy +1/2). Periodic boundary conditions 



13 




-1.0 

-1.2 ^ 

-1.4 £ 

I 

-1.6 1 — i 



-1.8 M 

o 



-2.0 



FIG. 9. Properties of the steady state as the strength £ of 
the OM driving is increased. The ratio between the number of 
particles in the symmetric and antisymmetric mode (dashed 
line, left axis) increases, while the trace of the separable ap- 
proximation p sep of Eq. (Bl I to the full density matrix (empty 



circles, right axis) approaches unity exponentially (solid line). 
We use ei = 0.2f, Ai = J = 0.2f , n = O.lf, and A = 2.0f. 
The state has been evolved to t max r = 10.0. The truncation 
of the Fock space is n max = 5 for the modes &e and n max = 1 
for the mode c, which is negigibly populated. 



are used in the Hamiltonian matrix. Now let us consider 
the transformation denned by the matrix 



M = 




(A2) 



such that MM T = 1 + 0[J 2 } and 



M II M 



-1/2 
1/2 1 

-1/2 

1/2 



+ 0[J 2 } . (A3) 



That is, we have presented a transformation M that, 
to first order in J/Slo, is orthogonal (i.e. implements a 
canonical basis transformation from independent bosonic 
modes to independent bosonic modes) and diagonalizes 
the hopping part of the Hamiltonian. In other words, the 
dominant part of the hopping is included into the defini- 
tion of the new basis states. The new basis states support 
the definition of the modes ag (with energy ~ —1/2) and 
eg (with energy ~ 1/2). The first correction includes hop- 
ping between next-nearest neighbors, that is between ag 
and &i±i or between ct and c?±\. 



Appendix B: Separability of the steady state 

In this appendix we discuss in more detail the effect 
of switching on the OM coupling in the two-membrane 
system introduced in Sec. Ill As we show in Fig. [9j 



the ratio betwen the population in the symmetric mode 
and the population in the antisymmetric mode increases 
monotonically with the strength of the driving. One 
has to keep in mind, though, that the strength cannot 
be increased indefinitely because the model that we use 
here adopts the rotating-wave approximation, which is 
not justified for arbitrary large drivings. The symmetric 
mode in the steady state is not only increasingly pop- 
ulated, but also increasingly disentangled from the re- 
maining part of the system, as the driving is increased. 
To demonstrate this, we produce an ansatz for density 
matrix in the steady state, in the form of a separable 
state p scp plus some remainder cr en t, and we show that 
the contribution of p sep to the total trace converges to 
unity exponentially. More precisely, we first diagonalize 
the steady-state density matrix as p = J^PaI^aX^aI, 
where p\ is a set of probability values that sums to 
unity and \^x) 1S a pure state. Each state |>Fa) can be 
expanded in the tensor basis of the modes as I^a) = 
Hn B n^n 3 *A ; n a n a n 3 l™s> l™a) K) , where n s , n a , n 3 arc the 
indices of the Fock states in the symmetric, antisymmet- 
ric, and c modes, respectively. We take n 3 = 0, corre- 
sponding to a negligibly populated excited state, and re- 
tain only the terms in the sum with the Fock occupation 



ni that maximizes J^ r . 
expression 



The resulting 



k)«|®|n 3 = 0}<n 3 = 0| 



(Bl) 



has the form of a separable state by construction, but 
non-unitary trace J2x S„ a ^l*A : n< A) ,n a ,ol 2 ' and tne ful1 
density matrix can be trivially decomposed as p = p sep + 
Cent- We see that, as the driving strength increases, the 
trace of p scp saturates the unitary total trace exponen- 
tially, meaning that the ansatz becomes numerically ex- 
act for strong driving. The separability of the steady 
state should be seen here as analogous to the purity of 
the steady state as stated in Eq. pi. Namely, the re- 
duced density matrix p s = Tr„ ain ^n a |(n 3 |p|n a )|7j3) of 
one mode is not a pure state both because of the en- 
tanglement between the modes, and because of the local 
coupling to the external reservoirs, which is effectively 
switched off in the number conserving models of driven- 
dissipative systems [BJ. In other words, the entanglement 
entropy, i.e. the von Neumann entropy of the reduced 
density matrix, is a good measure of entanglement for 
states that are pure, while we consider a system which is 
continuously subjected to pumping and losses. Separa- 
bility is, in this respect, the closest concept that we may 
apply here to characterize the fact that the combined 
process of driving and dissipation does not entangle the 
dark state of the system (the symmetric mode) to the 
other (auxiliary) modes. One may also expect that the 
steady state is actually a product state, i.e. a separable 
state where only one element in the sum over A is finite. 
We conjecture that this expectation is fulfilled in the limit 
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where the population of the states is very large, when the 
pumped antisymmetric state is substantially insensitive 
to the losses towards the symmetric state. Intuitively a 
larger intensity in the symmetric mode is correlated to 
a smaller intensity in the antisymmetric mode, because 
the coupling between the two modes is not negligible with 
respect to their internal dynamics. In the limit of large 
populations, however, the antisymmetric state can be ef- 
fectively seen as a bath that injects particles into the 
symmetric state, and hence the total system factors. 



Appendix C: Interpolation ansatz 



To validate the ansatz (37), one must make sure that 



the displaced thermal state is a solution to the equation 
£[p(a)] = for the steady state. Using the properties of 
the Glauber operators in the definition of p(a), we first 
compute the Liouvillian C^[a]p] = T>{—a)L[p{a)\D{a) 
in the displaced frame, which acts on the diagonal den- 
sity matrix p and where a is a parameter. The equation 
jCd [ce; p] — for the steady state is then treated within the 
coarse-grained approximation mentioned above. More 
precisely, in the equation of motion we neglect linear and 
cubic terms, which perturb the diagonal thermal state 
but do not substantially change the approximate steady- 
state manifold defined by the ansatz. The dominant part 
of the coherent fraction is contained in the displacement 
produced by the Glauber operators, while we neglect the 
quantum fluctuations of the coherent fraction induced by 
linear and cubic terms. Such quantum fluctuations are 
expected to be much faster than the slow dynamics of the 
order parameter a, and in this sense this approximation 
can be understood as a coarse-graining of the effective 



dynamics. The quartic terms vanish identically on the 
diagonal density matrix, while the remaining quadratic 
terms read 



C d [a;p] ~ R(n)[2apa) - {n,p}} 

+ f (n)[2a<pa- {aa f ,p}] 



(CI) 



This form consists of effective dissipation of width k(n) = 
jg)r / and incoherent pumping of 
|r, with V = ipr. The effect of 
the engineered dissipation oc Y is thus to renormalize the 
loss and pumping coefficients in the displaced frame. The 



n^-Y' 

strength f (n) = T+nf^Y', with Y' = ±fT 



right-hand side of Eq. (CI) vanishes if the thermal con- 
tribution to the density fulfills n — T(n)/[R(n) ~ T(n)], 
the solution of which gives 

(C2) 



K-T + \a\ 2 Y' ' 

The rescaled driving strength Y' includes model- 
dependent parameters so that the formula for n is simple. 
The large factor between Y' and Y is due to the com- 
plicated form of the mean-field Liouvillian, where many 
nonlocal terms add up. We verified that a different model 
Liouvillian, which describes onsite driving of photons be- 
tween bands in an array of equally spaced membranes, 
reduces as well to the present form, with different but 
equivalent functions k and T. 



Appendix D: Coefficients of the Gutzwiller 
approximation 

We derived the following convenient formula, which 
gives the mean-field decoupling for a Liouvillian in a non- 
diagonal Lindblad form, coupling up to M sites together 



T±*iWt, 



A I 



A I 



2 n 4t„ u n ^ 



A I 



A I 



n n 4 



A I 



(Dl) 



read 



with the definition of the local average (A)g — Tr^ Agpg . The coefficients of the mean-field Liouvillian in Eq. (35) 



2{2(a)[l + X + X 2 + 

-X(l + X)(a f a>] + 

+(4 + 7x + 4 X 2 )(a t aa)} 







2{2(at)[l + x + x 2 + 

-X(l + X)(a f a)] + 

+(4 + 7x + 4x 2 )(a t a t a)} 

4[{ato) +x(l + x)x 
x((at)(a} + (ata))] 

-2(x(« f ) 2 + 
+2(1 + X )(at 2 )) 
-4(l + 2 X + 2 X 2 )(at) 



-2( X (a> 2 + 
+2(1 + X)(a 2 » 

4(1 + X + X 2 )x 
x((ata> +1) 

2(2 + 3 X )(a) 



-4(l + 2 X + 2 X 2 )(a) 

2(2 + 3 X )(a t ) 
4(l + 2x + 2x 2 ) J 
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-2{(a)[(2 + 3 X ) + 

-2 X (l + x)(a t a)] + 

+ (4 + 7x + 4x 2 )(a t aa)} 







-2{<ot)[(2 + 3 X )+ 
-2 X (1 + X )(ata)] + 
-(4 + 7x + 4 X 2 )(atata)} 

4(1 + X + X 2 )(a f a) 

-2[x(at) 2 + 
+2(l + x)(a t2 >] 

2(2 + 3x)(a t ) 



-2(X(«> 2 + 
+2(1 + X )(a 2 » 

4[X(1 + X)(a f )(a) + 
+(1 + X + X 2 )x 

x«ata) + i)] 

-4(l + 2x + 2 X 2 )(a> 



2(2 + 3 X )(a) 

-4(l + 2x + 2 X 2 )(a t ) 
4(l + 2x + 2 X 2 ) / 
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